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O RADIATIVE TRANSFER IN PROTOPLANETARY DISKS 
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Ph. 

^^ ' Abstract. We present a new 3D continuum radiative transfer code, 

MCFOST, based on a Monte-Carlo method. The reliability and effi- 
ciency of the code is tested by comparison with five different radiative 
transfer codes previously tested by |Pascucci et al. (2004)| using a 2D 
^ \ disk configuration. When tested against the same disk configuration, 

no significant difference is found between the temperature and SED 
calculated with MCFOST and with the other codes. The computed 
values are well within the range of values computed by the other codes. 
The code-to-code differences are small, they rarely exceed 10% and are 

\^\ ' usually much smaller. 
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O |- 1 Introduction 
I 

^ ' Dust is present everywhere in astrophysics, from the interstellar medium to the 

"^ I atmospheres and close circumstellar environments of numerous classes of objects; 

d • from the lower mass planets and brown dwarfs, to massive stars. With the ad- 

^ [ vent of high angular resolution and high contrast imagers, the basic structural 

•"^ . properties (e.g., size, inclination, and surface brightness) of the circumstellar en- 

/\ ' vironments of the nearest and/or largest objects — disks and envelopes around 

i3 ■ young stars in nearby star forming regions and around more distant evolved stars 

■ ■ ■ ' — are now under close scrutiny. With this unprecedented wealth of high resolution 

data, from optical to radio, fine studies of the dust content become possible and 

potent radiative transfer (RT) codes are needed to fully exploit the data. 

At short wavelengths, dust grains efficiently absorb, scatter, and polarise the 
starlight while at longer wavelengths dust re-emits the absorbed radiation. How 
much radiation is scattered and absorbed is a function of both the geometry of the 
circumstellar environment and the properties of the dust. In turns, the amount 
of absorbed radiation sets the temperature of the dust (and gas) and defines the 
amount of radiation that is re-emitted at longer, thermal, wavelengths. 
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To get a reliable understanding of the structure and evolution of these "dusty" 
objects, be it the evolution of dust sizes, the temperature dependent chemistry, or 
simply the density profiles, it is highly desirable to model not only the integrated 
fluxes {i.e., the spectral energy distributions), but also the resolved brightness 
maps and polarisation profiles when available. This can only be done by solving 
the radiative transfer (hereafter RT) problem in media that can have large optical 
depths and/or complex geometries and compositions. Needless to say, analytical 
solutions do not always exist and sophisticated numerical methods must generally 
be used. One such versatile numerical method is the Monte Carlo method. The 
code we describe and test below is based on that RT scheme. 

The original version of the code, MCFOST, treated the polarised continuum 
RT problem by including scattering only, without a calculation of the temperature 
structure and thermal emission from the circumstellar environment Menard (1989) 



That code was used extensively to produce synthetic images of the scattered light 
from disks around young stars. Examples include the circumbinary ring of GG Tau 
Duchene et al. (2004)1 the large silhouette disk associated with IRAS 04158+2805 



(Menard et al. 2005), and an analysis of the circular polarisation in GSS 30 



Chrysostomou et al. (1997) In this paper we present a new improved version of 
MCFOST that includes passive heating of the circumstellar environment by stel- 
lar radiation and thermal re-emission from the dust, needed to prepare future 
observations with Spitzer, Hershel and Alma. 

To validate MCFOST, we perform a comparison with five other RT codes, in- 
cluding three using a Monte Carlo method that were benchmarked in a previous 



study Pascucci et al. (2004) (hereafter P04). To compare the reliability and ef- 
ficiency of each code, P04 used a well defined 2-dimensional disk configuration, 
with simple dust properties. Each code calculated the temperature structure in 
the disk, as well as the emergent SED. The results were compared quantitatively 
as a function of optical thickness and disk inclination. The agreement is gener- 
ally better than 10% for the most difficult cases, i.e., large optical depths, and 
much better in the optically thin cases. We performed identical calculations with 
MCFOST and present the results and comparisons below. 

In §2, we briefly describe MCFOST. In section §3, we recall the benchmark 
problem used to compare MCFOST with the 5 RT codes included in P04. Results 
are presented in §4, they include the temperature profile and SED for a range of 
disk optical thicknesses. 



2 Description of the numerical code 

MCFOST is a 3D continuum radiative transfer numerical code based on the Monte- 



Carlo method. It was originally developed by Menard (1989) to model the scat- 
tered light in dusty media (including linear and circular polarisations). In this 
paper, we present an extended version of the code that includes dust heating and 
continuum thermal re-emission. 
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2.1 Geometry of the computation grid 

The spatial grid is defined in cylindrical coordinates. It is well adapted to the 
geometry of circumstellar disks. We use Nr logarithmically spaced radial grid cells 
and Nz linearly spaced vertical grid cells. The size of the vertical cells follow the 
flaring of the disk, with a cutoff at 10 times the local scale height. For the purpose 
of the benchmark calculations presented in this paper, axisymmetry is assumed. 
However, MCFOST allows the density distribution to be defined arbitrarily in 
3D, with the limitation that within each cell, quantities are held constant. For 
example, MCFOST in a 3D version has already been used to the study the photo- 



polarimetric behaviour of AA Tau Pinte & Menard (2004) 



2.2 Position dependence of tlie dust distribution 

An explicit spatial dependence of the size distribution /(a, r) is implemented in 
MCFOST. The dust properties are therefore defined locally, i.e., each cell of the 
disk may contain its own independent dust population. This allows to model 
dust settling towards the midplane of the disk and/or variations of the chemical 
composition of the dust from the inner, hot regions to the outer, cold edge of the 
disk. 

2.3 Optical dust properties 

Dust particles can be arbitrary shaped and are assumed to be randomly oriented. 
In the case of homogeneous and spherical particles, the dust optical properties 
are computed with Mie theory. The optical properties in any cell of the disk are 
derived in accordance with the local size and composition distributions f{a,r), at 
the expense of a significant usage of memory space. The extinction and scattering 
opacities are given by 

«-*/--(A,f)- r""\a2Q,,t/sca(A,a)/(a,f)da (2.1) 

and the local MuUer matrix is 

S'(A,r)= / ""s'(A,a)/(a,r)da (2.2) 



where Qsca(A, a), Qcxt{^, a) and S{\, a) are respectively the scattering and extinc- 
tion cross sections and the Mueller matrix of a grain of size a at a wavelength 
A. 

2.4 The radiative transfer scheme in MCFOST 

The Monte Carlo method allows to follow individual monochromatic "photon pack- 
ets" that propagate through the circumstellar environment until they exit the 
computation grid. The propagation process is governed by scattering, absorption 
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and re-emission events that are controlled by the optical properties of the medium 
(opacity, albedo, scattering phase function, etc..) and by the temperature distri- 
bution. Upon leaving the circumstellar environment, photon packets are recorded 
into frequency dependent synthetic images and SED's can be constructed. 

To improve its computational efficiency, MCFOST uses different strategies to 
set the energy of a photon packet, corresponding to different samplings of the 
radiation field. The energy of a packet, hence the number of photons it contains, 
is chosen and optimised depending on the goal of the calculations. On the one 
hand, the convergence of the temperature distribution is optimized when all photon 
packets have the same energy, independently of their wavelengths. This procedure 
insures that more photon packets are emitted in the more luminous bins. On the 
other hand, the computation of SEDs is more efficient when it is the number of 
photon packets that is held constant for all wavelength bins. In that case it is 
the energy of the packets that is wavelength dependent. Thus, MCFOST is more 
efficient when it computes the temperature and SED with a two-step process: 

• Step 1 is the temperature determination. Photon packets are generated and 
calculated one by one. They are generated at the stellar photosphere and 
followed until they exit the computation grid. Upon scattering in the disk, 
the propagation vector of the packet is modified, but not its wavelength. 
Upon absorption however, packets are immediately re-emitted, in situ and 
isotropically, but at a different wavelength, calculated according to the tem- 
perature of the grid cell. For this re-emission process, the concept of immedi- 
ate reemission and the associated temperature correction method proposed 
by|Bjorkman fc Wood (2001) are used. They are presented as scheme 3 in 



P04. In this step, all photon packets have a similar energy and are ran- 
domly scattered/absorbed within the disk. The concept of mean intensity 



suggested by Lucy (1999) is further used to reduce noise in the temperature 



estimation for optically thin cases. Step 1 allows for a fast convergence of 
the temperature but is rather time-consuming when used to derive SEDs, 
especially in the low-energy, long wavelength regime. Therefore SEDs are 
calculated differently. 

• Step 2 computes the SED and/or images from the temperature distri- 
bution calculated in step 1. In step 2, the number of photon packets is 
held constant at all wavelengths and the enforced scattering concept of 



Cashwell & Everett (1959) is implemented to further speed up convergence. 
Step 2 maintains a similar noise level per wavelength bin and reduces conver- 
gence time for the SED by limiting the CPU time spent in high luminosity 
bins and focussing on low luminosity bins. 

2.4.1 The emission of photon packets 

In MCFOST, two radiation sources are considered, the star and the circumstellar 
environment. The stellar photosphere can be represented by (i) a sphere radiating 
uniformly or (ii) a limb-darkened sphere or (iii) a point-like source, whichever is 
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relevant for the problem under consideration. Hot and/or cold spots can be added 
the photosphere. In the following, wc consider a unique central point-like star that 
radiates like a black body. 

In step 1 of the previous subsection, all photon packets are emitted isotropically 
by the central star and have the same energy, independently of their wavelength. 
We define the luminosity of a packet, e, as: 

e = L,/N^,tepi (2.3) 

where L* is the bolometric luminosity of the star and N^stepi the number of packets 
generated. The wavelength of the packet is chosen according to a normalised 
probability density function proportional to B\{T^). 

In step 2, the SED is constructed using a constant number N^step2 of packets 
for all the wavelength bins. The packets are emitted by the star and by the disk. 
The energy of the N^step2 packets emitted at a given wavelength A is determined 
by the total energy that the star and the disk radiates at this wavelength. The 
corresponding packet luminosity is : 

X.(A) + E, L,{\) 

CA = ^ (2.4) 

where i is the index of the cell and i*(A) and Li{X) are the luminosities of the 
star and disk cell at a given wavelength. They are given by: 

L4\)=iTr^RlBxin) , (2.5) 

and L,{X) = 47r nf^X) Bx{T,)m^ . (2.6) 

The wavelength is set in a deterministic way and packets are randomly emitted 
from the star and the disk, by cell i, with the respective probabilities: 

P* L,(A) + E. i.(A) ' ^'■'> 

-'^^'^ x.(a)+'e.l.(a) - ^'-'^ 

In each disk cell we assume that the density, temperature and opacities are 
constant. To determine the position of emission of a photon within a given disk 
cell, in cylindrical geometry, the following relations are used: 



mill ' ^v max ' min-' 



Mrl.. - rlJ , (2.9) 

\z\ — ^min + ^(^max ^ 2;fnin) j (2-10) 

and (p = 2TrA (2.11) 

where A denotes three different random numbers drawn from a uniform distri- 
bution in the interval ]0,1] and where r^in, Tmax: -Zmin: 2^max are the radial and 
vertical boundaries of the cell. The dust thermal emission is assumed isotropic. 
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2.4.2 Distance between interactions 

It is natural within a Monte Carlo scheme to estimate the distance a photon packet 
"travels" between two interactions by means of optical depth (related naturally to 
the density) rather than by physical, linear distance. From a site of interaction, 
the optical depth to the next site is randomly chosen from: 



T\ 



In A 



(2.12) 



with A G ]0, 1]. The distance I is computed by integrating the infinitesimal optical 
depth K°'^*(A, r)p{r)As until the following equality is verified : 



T\ 



K'='<^*(A,f)p(f)ds 



(2.13) 



Once the position of interaction r is determined, the probability that the in- 
teraction is a scattering event rather than an absorption event is estimated with 
the local (i.e., cell) albedo : 



Psca = ^(A, r) = 



J 7ra^Qsca(A, a)f{a, r)da 
J 7ra2Qoxt(A, a)f{a, r}da 



(2.14) 



In step 1, the photon packets are scattered or absorbed following this probability 
whereas in step 2, they are always scattered, no absorption is allowed to occur 
explicitely. Instead, scattering is enforced at each interaction site but the Stokes 
vector is weighted by the probability of scattering, Psca- This allows all the photons 
to exit the disk and to contribute to the SED, although with a reduced weight. 



2.4.3 Scattering 

MCFOST includes a complete treatment of polarisation, including linear and cir- 
cular polarisations. In the Stokes formalism used, the state of a light packet is 
described by its Stokes vector, with / the intensity, Q and U describing the linear 
polarisation, and V the circular polarisation. The interaction of a photon with a 
dust particle is described by a 4 x 4 matrix, the Mueller matrix, S. In the simpli- 
fying case where dust particles are homogeneous and spherical (Mie Theory), the 
matrix becomes block-diagonal, only 4 different elements are non-zero, and only 
three of these are independent. 



Q 

u 



+1 



/ 



V 



5ii 

'S'l2 






5*12 

5*11 
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\ 



( M 
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u 

\y J 



(2.15) 



The direction of scattering is defined by two angles 9 and (j) in spherical coordinates. 
Each individual element Sij of the matrix is a function of the two scattering angles 
and of the wavelength, Sij — Sij{0, <p, A) 
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The scattering angle 9 is randomly chosen from the pretabulated scattering 
phase function : 

^^i0)^Sn{0) (2.16) 

where 6*11 is the first element of the Mueller matrix. 

For unpolarised incoming light packet, i.e., Q,U = 0, the distribution of az- 
imuthal angle is isotropic. For light packets with a non-zero linear polarisation 
P = ^JCp^^^^lP- 1 1 , the azimuthal angle is defined relative to its direction of polar- 
isation and determined by means of the following repartition function : 

p.,. 1 /, S,m~S,M sin(20) \ 



where Q was previously chosen from equation 12. l()l (see [Sole (1989) I 



The Mueller matrix used in equation 12.1 51 [TT7I and !?. 161 can be defined in two 

different ways : one Mueller matri xper grain size or one mean Mueller matrix per 

grid cell. In the latter case, the Mueller matrix is defined by eouation l2.2l In the 

former case, the grain size must be chosen explicitely for each event following the 

probability law : 

7ra2Q,^a^(a)/(a, f)da 

p[a)aa = —^ ■r~r~r ^ — • (2.18) 

JaZT T^a Qsca(a)/(a, r)da 

The two method are strictly equivalent. They can be used alternatively to optimise 
either the memory space or the computation time required. 

2.4.4 Absorption and radiative equilibrium 

The temperature of the dust particles is determined by assuming radiative equi- 
librium and that the dust opacities do not depend on temperature. 

Each time a packet 7 of a given wavelength A travels through a cell, we compute 
A^, the length travelled by the packet in the cell and the mean intensity J\ in 



the cell is derived following Lucy (1999) 



Jx = ^—YeAL^^—Y — AL (2.19) 

1 1 ' 

where Vt is the volume of the cell i. 

The dust thermal balance should take into account the thermal coupling be- 
tween gas and dust. In high density regions, close to the disk midplane, the 
coupling is very strong and the gas temperature should be close to the dust tem- 
perature. In the surface layers of the disk, on the other hand, density becomes 
very low and gas-dust thermal exchanges should be reduced. 

The two extreme assumptions are implemented in MCFOST : either the gas- 
dust thermal exchange is perfectly efficient and gas and dust are in local thermo- 
dynamic equilibrium (LTE), either there is no thermal coupling between gas and 
dust. 
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If we assume local thermodynamic equilibrium, all dust grains have the same 
temperature which is equal to the gas temperature. We can define this temperature 
as the temperature of the cell, given by the radiative equilibrium equation : 

/•oo 1*00 

47r / nf'\X)Bx{T,)d\ = ATT Kf'{X)JxdX (2.20) 

which we write : 

j^ nf%\)B,{T,)dX =—l—^Y.^f^iX)Al, . (2.21) 

If there is no thermal coupling (other than radiative) between gas and dust, and 
then between the different dust grain sizes, each grain size has its own temperature 
and the radiative must be, in this case, written independently for each size : 

/•OC 



4n Kf''{\,a)BxiT,{a))d\^ATT Kf%X,a)JxdX (2.22) 

Jo Jo 

where Kf^^{X,a) and Ti{a) are the opacity and temperature of the dust grains of 
size a in the cell i. 

The calculation of J^ Kf^^{X)Bx{Ti)dX = aT'^Kp/n {kp is the Planck mean 
opacity) is very time consuming and we pretabulate these values at Nt = 1000 
logarithmically spaced temperatures ranging from IK to 1500K. The temperature 
Ti of each cell is obtained by interpolation between the tabulated temperatures. 



Following Bjorkman & Wood (2001) we use the concept of immediate reemis- 



sion. When a packet is absorbed, it is immediately reemitted and its wavelength 
is chosen taking into account the temperature correction following the probability 
distribution : 

pxdX OC Kf%X) f^\ dX. (2.23) 

3 Benchmark problem 

In order to test MCFOST, we calculated the same cases as described in P04 and 
compared the results. The geometry tested involves a central point-like source 
radiating as a T=5 800K blackbody encircled by a disk of well defined geometry 
and dust content. The disk extends from 1 AU to 1000 AU. It includes spherical 
dust grains made of astronomical silicate. Grains have a radius of 0.12 /im and 
a density of 3.6 g.cni"^. The optical data are taken from Draine & Lee (1984) 



The disk is flared with a Gaussian vertical profile p{r, z) = po{r) exp{—z'^/2 h{rY) 
Power-law distributions are used for the surface density S(r) = Eg (r/ro)" and the 
scale height h{r) — Hq (r/ro)^ where r is the radial coordinate in the equatorial 
plane, ho the scale height at the radius tq. These assumptions lead to a general 
expression for the density at any point in the disk : 

\ a-l3 / , ■i^ 

r \ 1 



pir,z)=po[-J cxpl--(^— jl. (3.1) 
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In this benchmark, /3 — 1.125, a ~ (3 — —1.0 are used. The reference radius 
is tq = 500 AU, the disk height is Hq = 99.74 AU. This scale height is the same 
as described in P04. A factor \f^lpK appears because of the different definition of 
the Gaussian profile in MCFOSTi. 

The scattering is set to be isotropic and polarisation is not calculated, i.e, all 
the scattering information is contained in the scattering cross section alone, Qscai 
acting on the I Stokes parameter only. 

Results are shown in the section below for 4 different disk equatorial optical 
thicknesses, ry : 0.1, 1, 10, and 100. For each optical thickness case, results for 
three different disk inclinations are calculated, 12.5°, 42.5° and 77.5°. 



4 Results 

4.1 Computational considerations 

For the test cases considered here, the number of grid cells is set to Nr = 50 
and Nz = 20 in the radial and vertical directions. 10^ photon packets are used 
to calculate the temperature distribution (step 1) and 2 10® photon packets per 
wavelength are used for the generation of the SED (step 2). The total runtime is 
13 minutes for the most optically thin case and 20 minutes for the most optically 
thick case on a bi-processor (Intel Xeon) computer running at a clockrate of 2.4 
GHz. The runtime memory space needed is 10 MBytes. 

For the high resolution models (see sec. 14.51 below), with Nr — 500 and N^ = 
200, the runtime for the most optically thick case is « 7 hours and the memory 
occupation is 450 MBytes. 



4.2 The temperature profiles 

The radial and vertical temperature profiles for the most optically thick case are 
shown in the top panels of Fig. ^ and El For clarity the results of MCFOST are 
shifted by 200 K in Fig. ^and by 40 K in Fig. [2 The lower panels of both figures 
show the difference, given in percents, taking the code RADICAL as a reference. 
The thick sohd lines traces the difference of MCFOST with RADICAL. The radial 
temperature of MCFOST does not differ by more than 5% from all the other codes, 
except from MC3D and RADICAL close to the inner radius of the disk, and from 
STEINRAY at large radius, but the maximum difference remains below 15%. In 
the vertical direction, MCFOST always agrees better than 2.5% with RADICAL, 
MCTRANSF and RADMC. The agreement with MC3D and STEINRAY is always 
better than 2.5% at high altitudes {9 > 20°). Closer to the midplane deviations 
are larger but do not exceed 4%. 



^See the definition of /2(r) in eq. 4 of P04. 
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Fig. 1. Plots of the radial temperature (upper panel) and the difference (lower panel) 
for the most optically thick case, rv ~ 100 and using RADICAL as the reference. In 
both panels, the results of MCFOST are represented by the thick solid line. Thin solid 
lines are the results from MC3D, dot-dashed lines from MCTRANSF, dashed-dot-dot-dot 
from RADICAL, dotted hues from RADMC and dashed hues from STEINRAY. In the 
upper panel, and all curves being very similar, MCFOST has been shifted by 200K for 
clarity. 



4.3 The spectral energy distributions 

Shown in Fig. Oare the calculated SEDs for two disk inclinations and four different 
optical thickness, Ty. The figure directly compares with Fig. 7 of P04. We plot 
XFx in (W.m~^) where Fx is the flux density at a distance equal to the stellar 
radius, for a naked star (without disk) Fx — ttBx (triangles in Fig. ^. MCFOST 
reproduces the correct slope at long wavelengths, i.e., XFx oc A~^, expected for an 
optically thin medium containing small particles (with k oc A^^), emitting at long 
wavelength (Ba oc A"'*). 

In Fig. |5l the difference of the SEDs calculated by the various RT codes are 
presented. The overall agreement is better at lower inclinations, a consequence of 
the lower optical depth along the observer's line of sight. Indeed, for models with 
Ty = 0.1 and Tv = 1 at all inclinations and for models with ry = 10 and Ty = 100 
at inclinations of 12.5° and 42.5°, deviations between MCFOST and all the other 
codes do not exceed 10%. For the most optically thick cases, i.e., large inclination 
and/or large optical thickness, the largest differences are observed below 1/im, 
where scattering dominates, and around lOyitm, in the silicate band. 
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r = 2 AU 




Fig. 2. Plots of the vertical temperature (upper panel) and the difference (lower panel) 
among the codes using RADICAL as the reference, for the most optically thick case, 
TV = 100 and for a distance r in the midplane equal to 2 AU from the central star. In 
the upper panel, the results of MCFOST are shifted by 40 K for clarity. The line types 
are the same as in Fig. 



In the optical range, MCFOST produces slightly larger fluxes with respect to 
other codes but the large dispersion observed in all codes does reflect lower signal 
to noise ratio due to lower emerging fluxes. Near 10/im, MCFOST shows the same 
trend as MCTRANSF but flnds larger fluxes than RADICAL or MC3D. 

To accelerate convergence in MCFOST, photon packets are recorded in incli- 
nation bins that cover a range of angles, instead of a single angle as is often done 
elsewhere. It is important to note that by doing so, the uncertainty associated 
with the models easily matches the observational errors. Observationally, inclina- 
tion angles of disks are rarely known to better than 5-10 (1-2) degrees when seen 
pole-on (edge-on). In MCFOST, the inclination bins intercept equal solid angles 
from pole-on to edge-on, i.e., they cover equal intervals in cosj. This interval 
can be set to match the quality of the data. In the calculations presented in this 
papers, 21 bins are used from pole-on to edge-on. In all flgures, the inclinations 
listed, i.e., 12.5°, 42.5°, 77.5°, are the median of the bins covering an interval of 
0.05 in cosi. These bins have inclinations ranging from 0° to 17.75°, from 40.37° 
to 44.42°, and from 76.22° to 79.02°, respectively. 

For comparison, the axis ratio of a circular disk estimated with a measurement 
error equal to half the bin width used in this paper, i.e., 0.025, would yield an 
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1000 



Fig. 3. SED for a disk inclination i — 12.5°. Each curve traces a model SED calculated 
by MCFOST. The midplane optical depth is given in parenthesis above each curve. The 
solid lines show results for the most optically thin disk, t^, — 0.1, dotted lines for a 
disk having t^, — 1, dot-dashed lines for a disk with r^ = 10, and dashed lines for the 
most optically thick model, r^ — 100. The diamonds trace the black-body emission from 
the naked star. The slope of the SED at long wavelengths depends only on the dust 
properties and is plotted in each panel with a solid line, Fa oc A~^. 



error on i of ~ 10° for the pole-on case, and ~ 1.5° in the edge-on case. Indeed, 
for a perfectly circular pole-on disk, measuring an axis ratio of 0.975 instead of 
1.0 would mimic an inclination of i = 12° instead of i = 0°. In the edge-on 
case, measuring an axis ratio of 0.1 instead of 0.125 would induce a difference of 
1.4° in the determination of the inclination. The results presented here, generally 
obtained using 21 inclination bins, are therefore accurate enough and not limited 
by our choice of inclination bin number. 



4.4 Impact of inclination angle sampling 

To verify further that an inclination range (as used) does not introduce a noticeable 
bias in our comparisons, we show in Fig. |51 the same calculations for inclination 
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10 100 

Fig. 4. Same as Fig. I^for an inclination angle of i = 77.5° 



bins 20 times narrower. In Fig. the bins run from 11.76° to 12.89°, 42.31° to 
42.67°and from 77.38° to 77.63°. The results are similar to those of Fig.El bottom 
right panel. 

4.5 Impact of grid sampling 

In this section, we test the influence of the grid sampling on the resulting SEDs. 
In Fig. ^1 the difference between the previous calculations with Ty = 100 and the 
same calculations performed with a grid sampling 10 times finer, both in radial 
and vertical direction are shown. The number of the photons is the same for 
the two simulations. The high resolution calculations were done with the same 
number of photons as previously, the corresponding runtime is ~ 7 hours. The 
comparison shows that the spatial resolution has very little influence on the SED, 
with differences not exceeding 2.5% even in the worst case. 



4.6 Impact of vertical cutoff 

MCFOST is the only code which uses a cylindrical grid and introduces a verti- 
cal cutoff for the density, all other code have adopted spherical geometries and 
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Fig. 5. Plots of the difference between the model SED, given in percents, using the code 
RADICAL as the reference. Results from all codes are shown: the five codes included 
in the study of P04 and MCFOST (this study). Results are shown as a function of 
optical thickness of the disk midplane and as a function of inclination angle. Results are 
shown for the most optically thin case r^ — 0.1. Three different inclination angles are 
considered: i = 12.5° (upper), i = 42.5° (middle), and i — 77.5° (lower). All codes are 
compared to RADICAL. The thick solid line gives the differences for MCFOST, the thin 
solid line for MC3D, the dot-dashed line for MCTRANSF, the dotted line for RADMC, 
and the dashed lines for STEINRAY. 



do not need to set up a vertical cutoff. The cutoff is generally chosen at 10 
times the local scale height of the disk. The corresponding density is equal to 
exp(— 10^/2) w 210^^^ times the density in the midplane of the disk and is so 
low that no absorption or scattering actully occurs in this region. To verify the 
influence of the cutoff on the emerging spectra, we run the previous models with 
a cutoff of 20 times the scale height rather than 10, keeping the same vertical res- 
olution and then using twice more cells. The difference between the models with 
a cutoff at 10 and 20 times the scale height are shown in Fig. ^J They always 
stay below 3 % even at the highest optical depth. 
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Fig. 6. Same as figurc|^but for r^ = 1. 



5 Discussion and conclusions 

In the sections above we presented a new continuum 3D radiative transfer code, 
MCFOST. The efficiency and rehabihty of MCFOST was tested by considering a 
simple and weh-defined problem, a dusty disk as described in P04, that was used as 
a benchmark to compare MCFOST with 5 other RT codes. MCFOST was shown to 
calculate temperature distributions and SEDs that are in excellent agreement with 
previous results from the other codes. The results of MCFOST for all the test cases 
are available at the web site : http: //www-laog. obs .ujf-grenoble . f r/~cpinte/ mcf ost/| 

Because MCFOST reproduces well the results presented in P04, their conclu- 
sions also hold for MCFOST. In particular, we find that the far-infrared part of 
the SEDs is unaffected by the viewing angle, see Fig. 13 On the other hand, 
the spectrum at shorter wavelengths is largely modified, with the stellar spectrum 
dominating for pole-on views and becoming progressively more attenuated with 
increasing viewing angle. For cases with large attenuation of the starlight (high 
inclinations and large optical depths) the contribution from scattering becomes 
important because the dust grains have a high albedo. Also, the spectroscopic 
signatures of silicates at 10/im and 20/im are visible, as expected. The overall 
temperature profiles in the disk and shape of the emerging SEDs are well repro- 
duced by MCFOST. Differences in the radial temperatures computed by MCFOST 
do not exceed 4% with respect to the other 5 codes in the optically thin case. For 
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Fig. 7. Same as figurel^but for r^ = 10. 
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the most optically thick case, differences do not exceed 15%. Regarding SEDs, 
MCFOST is similar to all other codes to better than 10% for low optical thick- 
ness. Devitations as large as 15% are however observed in the most defavorable 
case, i.e., high tilt and optical thickness, except with STEINRAY which deviates 
from all other codes, between 10 and SO/im (see P04 for more detailed description 
of these deviations). 

In this paper, disks with fairly low equatorial V-band optical depths were 
calculated, ry = 100. To reproduce observed edge-on disks, MCFOST has been 
able so far to handle more compact, equatorially concentrated disks with ry '^ 10^, 
illustrating the ability of MCFOST to model much more optically thick disk than 
the one used for the benchmark simulation 

Various tests were further performed to show the independence of the MCFOST 
results with respect to grid sampling, inclination sampling, and position of the 
vertical density cut-off in the disk. Results are not noticeably affected. 
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Fig. 9. Plots of the difference of the results MCFOST and other codes for the most 
optically thick case with high resolution on inclination bins. The figure presents the same 
calculations as the bottom-right panel of Fig. |^but with 20 times narrower inclination 
bins. The lines type are the same as in Fig. |3 
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Fig. 10. Spatial resolution tests. Difference percentage between the emerging SED of 
the model with a grid of (500 x 200) points (radial x vertical) and the reference model 
with a grid of (50 x 20) points. The number of "photons packets" is the same for the two 
models. 
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Fig. 11. Effect of the vertical cutoff. We plot the difference, given in percents, of the 
model with a vertical cutoff at 20 times the scale height relative to the reference model 
plotted in|3 (cutoff at 10 times the scale height) for the most optical thick case. The two 
models have the same vertical resolution. 



